Stochastic Mean-Field Theory: Method and Apphcation to the Disordered 
Bose-Hubbard Model at Finite Temperature and Speckle Disorder 
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We discuss the stochastic mean-field theory (SMFT) method which is a new approach for describ- 
ing disordered Bose systems in the thermodynamic limit including localization and dimensional 
effects. We explicate the method in detail and apply it to the disordered Bose-Hubbard model at 
finite temperature, with on-site box disorder, as well as experimentally relevant unbounded speckle 
disorder. We find that disorder-induced condensation and reentrant behavior at constant filling are 
only possible at low temperatures, beyond the reach of current experimentsi^^. Including off-diagonal 
hopping disorder as well, we investigate its effect on the phase diagram in addition to pure on-site 
disorder. To make contact to present experiments on a quantitative level, we also combine SMFT 
with an LDA approach and obtain the condensate fraction in the presence of an external trapping 
potential. 

PACS numbers: 67.85.Hj, 03.75.Hh, 71.55.Jv 



I. INTRODUCTION 



The interplay between disorder and interactions in 
Bose systems has been a vital field of research in con- 
densed matter both in theory and experiment. The line of 
investigation was mainly initiated by the seminal work of 
Fisher et alji, who first provided a detailed study of local- 
ization of interacting bosons in a random potential, which 
led to the notion of the supcrfluid-insulator transition 
and the Bose glass (BG). While disorder effects in Fermi 
systems are relevant to a broad range of experimentally 
accessible scenarios like correlated electron systems, the 
status was less diverse for Bose systems for a considerable 
time period, as superfluid '''He situated in random pores 
of Vycor had been the predominant setup which could be 
studies with sufficient precision^. This changed dramati- 
cally when the realization of the supcrfluid-Mott insula- 
tor transition of ultracold bosonic atoms in an optical lat- 
tice opened up a new field of investigational^. In particu- 
lar, optical lattices provide a relatively pure and tunable 
simulation of effective models used to describe solid state 
systems^, where effects like disorder can also be realized 
in a controlled manner. While several alternative real- 
izations of disorder in optical lattices, such as multichro- 
matic lattices with non-commensurate wavelengths^^— 
or multi-species gases with strongly differing tunneling 
ratesiiii^ have been proposed recently, speckle laser pat- 
terns are probably by now one of the most efficient meth- 
ods to establish disorder in cold atomsH^— Therein, it 
is possible to combine the speckle beam with the remain- 
ing apparatus of the optical lattice to simulate disordered 
lattice systems with a high tuning accuracy and without 
other side effects. 

A variety of theoretical approaches^ has by now been 
applied to the disordered Bose Hubbard model (BHM), 
first introduced for ultracold atoms by Jaksch et ali^, 



which is described by the Hamiltonian 

Hbh = -jJ2iblb^+h.C.) 

{id) 

+ l-i)ni + j'^ni{n,-l), (1) 

i i 

where (&|) annihilates (creates) a particle in the low- 
est band Wannicr state at site i, ni ~ b\b^ is the local 
particle number operator, J denotes the nearest neighbor 
hopping energy in the lowest band, and ^ is the chemi- 
cal potential, ti is an on-site energy shift, which in our 
case is a spatially uncorrelated random variable drawn 
from a distribution p(e) and U is the on-site repulsive 
interaction. The subscript indicates the sum over 

all neighboring pairs of sites. Unless stated otherwise, we 
use the unit J7 = 1. 

Several quantum phases can exist within this model, 
such as the Mott insulator (MI), the Bose glass, the con- 
densed phase, commonly referred to as the superfluid 
(SF), as well as the normal phase at finite temperature. 
The transitions between these phases, which constitute 
some of the first experimentally feasible quantum phase 
transitions in bosonic systems, have attracted much at- 
tention. Numerically, a powerful approach is Quantum 
Monte Carloi^-^ (QMC). While MI and SF phase can be 
characterized efhcientlyi^, the BG phase and the vicini- 
ties of the transition lines are significantly more com- 
plicated to be adequately described. The main reason 
is that for finite size calculations in general, it is prob- 
lematic to capture the correct description of the phase 
borders, which are essentially dominated by rare events, 
which is also a problem for QMC methods. In most cases, 
exact diagonalization studies are simply inadequate due 
to limited size and number of particles, which often ob- 
scures essential physics (however, there may appear as- 
pects that can indeed be suitably captured by small clus- 
ters'^). On the other hand, with a similar range of 



treatable system sizes as QMC, density matrix renor- 
malization group (DMRG) is an efficient complementary 
method^I. However, DMRG is currently only applicable 
in one spatial dimension and thus does not allow for a 
description of effects in higher dimensional lattices. An- 
alytically, renormalization group analysia^Sri^, slave bo- 
son theory2i^ the strong coupling approach^, and vari- 
ous different kinds of mean field theory descendants have 
been applied to the BHM with and without disorder—"—. 
Arithmetically averaged mean-field theories, on the other 
hand, arc incapable of resolving the BG phase at all for 
T = 0^, but also for T > impose an unphysically 
strong phase coherence, leading to an overestimation of 
the SF phase. Other methods like the random phase ap- 
proximation (RPA) which has successfully been applied 
to the system without disorde r'^^-'^^ again suffer from fi- 
nite size effects for the BHM with disorder, as the absence 
of translational symmetry constrains its applicability to 
much smaller system sizes. 

To circumvent this type of problems in describing the 
different phases of the disordered BHM, we use the SMFT 
which has been previously introduced and applied to the 
disordered BHM at zero temperature^. There, it was 
found that the method efficiently describes localization, 
is both valid in high dimensions and in the thermody- 
namic limit, capturing rare events with their respective 
statistical weight and includes dimensional effects. In 
particular, it was found that at fixed /x, there exists a crit- 
ical hopping strength, above which the system remains 
supcrfluid for arbitrarily strong disorder. 

In this article, we present the SMFT in detail and in- 
vestigate how the results found for the disordered BHM 
at zero temperature are modified at finite temperature. 
In addition to prototypical box disorder, we consider ex- 
ponential speckle disorder to better simulate systems re- 
alized in current experiments. Here, the results are qual- 
itatively different in the sense that for any finite disorder 
strength the MI gives way to the BG. 

The paper is organized as follows: In Sec. [Hi the SMFT 
is explained in detail. First the general scope is outlined, 
followed by the definition of quantities computed within 
SMFT, such as compressibility, local Greens functions 
and condensate fraction. In Sec. IIIIl the essential results 
for the disordered BHM at T = are given, followed 
by the extension of the SMFT calculations to T 7^ 0. 
In Sec. IIVI we extend the results presented for box dis- 
order—, including finite temperature effects. Alterna- 
tive types of disorder, in particular disorder induced by 
speckle lasers arc discussed for the BHM in Sec.|Vl As an- 
other possible source of disorder, we discuss the effect of 
kinetic (hopping) disorder in the BHM in Scc. lVBl How- 
ever, we find no sensible dependence of the system on this 
parameter. The applicability of these results to current 
experiments essentially relies on the estimate of experi- 
mental temperature, which is discussed on in Sec. IVI Al 
We find that the experimentally realized temperatures 
are still far above the regime for which we resolve the pre- 
viously stated interesting phenomena, such as disorder- 



induced condensation and reentrant behavior. LDA -|- 
SMFT calculations are discussed in Sec. IVI Bl to provide 
a closer connection to experimentally measurable quan- 
tities. In Sec. IVIIl we conclude that the SMFT is an 
efficient theoretical approach to disordered Bose systems 
and promises an adequate description of ongoing experi- 
ments. 



II. METHOD 

As pointed out previouslyii^i, performing a self- 
consistent disorder average over all local on-site ener- 
gies is not sufficient to generally describe the insulat- 
ing Bose-glass phase. From spatially resolved bosonic 
Gutzwiller calculations, it becomes apparent that this 
method overestimates long range correlations, predicting 
the formation of a global condensate into a single parti- 
cle orbital which consists of the superposition of a large 
(extensive) number of distinct localized single particle 
states. In the true Bose glass phase, off-diagonal long 
range order does not prevail and a large (although not 
extensive) number of particles may occupy each of these 
localized modes independently, leading to a condensate 
fraction /c = in the thermodynamic limit. However, 
by imposing phase rigidity and averaging over all mean- 
field parameters ipi = (bi) with the same complex phase, 
the spatially resolved (as well as the arithmetically aver- 
aged) Gutzwiller theory leads to a finite average mean- 
field parameter (MFP) and a finite condensate fraction 

fc " (b) I {b^b) in the expected BG regime. 

In contrast to the approaches mentioned above, the 
SMFT is constructed as a single-site theory in the ther- 
modynamic limit, which effectively describes fluctuations 
in the MFFs using a probability density function (PDF) 
P{ip). Set up in this fashion, the SMFT method is not re- 
stricted to incorporating disorder fluctuations only, but 
may also be a powerful approach to treat fluctuations 
of different origins, such as of thermal or quantum type 
in a unified framework. The concept of approximating 
quantum mechanical operators by random variables has 
been applied previously in a variety of physical scenar- 
ios^"—. In this paper, we will focus on disorder-induced 
fluctuations at flnite temperature. In this section, we will 
discuss the construction of SMFT in detail, concentrat- 
ing on the case of including on-site disorder. Extensions 
of including disorder in the hopping parameter J or ther- 
mal fluctuations explicitly are discussed in Sec. IV Bl and 
App. [F] respectively. 

The central quantity, which effectively describes a dis- 
ordered bosonic lattice gas in the thermodynamic limit 
is the probability distribution P{ip). It is assumed to 
be equal for all sites and in particular independent of 
the nearest neighbors' on-site energies. The validity of 
this assumption has been checked using spatially resolved 
mcan-fleld theory, which is known to become exact in 
the non- interacting limit, while retaining interaction ef- 
fects beyond Gross-Pitaevskii theory for flnite U/J. It 
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is found that it yields correlation coefficients below 0.05 
in all regimes considered^, justifying the above assump- 
tion. The aim of SMFT is to find a self-consistent solu- 
tion for P{ip), which is restricted and uniquely specified 
by self-consistency equations and minimization of energy. 

Consider a cluster of lattice sites composed of a central 
site i and Z (the coordination number) nearest neighbors, 
where the corresponding part of the Hamiltonian ([T} con- 
tains the operators on site i within the bosonic Gutzwillcr 
approximationi*^^— : 



+ {e,-^,)blb, + ^blblbA- 



(2) 



Within this approximation, the site i is coupled to (and 
the sum extends over) the nearest neighbors via the scalar 
MFPs ijjj = (^i)j which arc random variables within 
SMFT. Due to global particle number conservation, the 
MFPs can all be chosen to be real and positive, as any 
variation in complex phase corresponds to a boost in lo- 
cal kinetic energy. The expectation value (bj) is to be 
evaluated in the local ground state for T = or taking 
the thermal trace for T > 0. Thermal fluctuations can, 
however, also be incorporated on an explicit stochastic 
level in the distribution P{tp) within SMFT, as discussed 
in App. If] Inspecting the Hamiltonian ([2]), it does not 
depend on the ipj's individually, but only on the scaled 



n.n.j 



(3) 



Since the tpj^s (and in the hopping disorder Sec. IV Bl also 
J) are random variables, the newly defined quantity r]i is 
also a random variable obeying some distribution Q{ri). 
As discussed above, the random variables ^pj are assumed 
independent within SMFT, allowing us to express the 
distribution Q{ri) by a Z-fold scaled convolution 



Qiv) 



?n — 1 



(4) 



Making use of the convolution theorem, this can be re- 
duced to two one-dimensional Fourier transforms by in- 
troducing the characteristic function 

^(t) = j d^P(^)e'*'^, (5) 

in terms of which the function Q can be expressed as 

= ^ JdtMt)f e-^'^/-'. (6) 

This can be calculated efficiently for arbitrary coordi- 
nation numbers Z using the FFT algorithm. The numer- 
ical procedure is discussed in App. 



The self-consistency condition can now be formulated 
in the following way: If the on-site energy e is randomly 
drawn from p(e) and ^pj is randomly drawn from the 
self-consistently determined distribution P{ip) for each 
of the Z nearest neighbors, this defines the single site 
Hamiltonian which can be diagonalizcd providing 
the ground state |g.s.(e, ry)). From there, a new MFP 
(g.s.(e, ?7)| 6 |g.s.(e, 77)) can be calculated, with the self- 
consistency requiring the distribution of this new random 
variable to be exactly the distribution P{ip) we initially 
assumed for the neighboring ■0j's. 

To cast this condition into functional form for a prob- 
ability distribution in the thermodynamic limit (i.e. for 
an infinitely large system at fixed density) , we first define 
the function 



5(m - e, V) = (g-S-(M -'^,v)\b |g.s.(Ai - e, v)) , (7) 

where |g.s.(/i — e, 77)) is the ground state of H,-^"'' (/i — e, rj) 
given in ([2|). 

It is useful to introduce the conditional PDF, which 
is a function of rj and "0, giving the probability density 
for a specified value ip if the value of ij is fixed and e 
is distributed according to p(e). This can be obtained 
by using the transformation property of a PDF under a 
variable transform and takes the form 



d 

~ dip 



de p{e) Q {ip - g{fi ~ €,7])) . 



(8) 



This function does not obey any self-consistency condi- 
tion and can be evaluated directly, which is the first step 
in finding P{ip)- Remembering that the relation between 
Q{r]) and P{ip) is given by so that we can express 
the self-consistency condition as 



P{ij)= / d77Q(7?)P(^|77). 



(9) 



The right hand side can be understood as follows: for 
every fixed rj we have a PDF P{ip\ri), which yields a 
contribution with the respective weight Q{r]), leading to 
a marginal distribution for the considered site. If this 
agrees with the initially assumed distribution P{'ip), a 
self consistent solution has been found. 

Once this is determined, expectation values of local, 
self-averaging operators A can be directly be determined 
by 



(i)=Tr(e(/3)) 



where 



dep{e) / dijQir]) 



g-,3«<"^'(M-e,„) 

Trfe-'3«^"^'(M-^,'7)' 



(10) 



(11) 
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is an effective disorder averaged density operator, incor- 
porating thermal, on-site energy and MFP fluctuations, 
depending explicitly on p{e) and the self-consistently de- 
termined Q{vi). 



A. Numerical solution 

To solve the SMFT equations numerically, we iterate 
the self-consistency equations on a discrctized grid for 
■0, consisting of a superposition of a variable number of 
equidistantly spaced grids, as explained in App. [X] For 
every fixed set of physical parameters, we first numeri- 
cally determine the conditional cumulative density dis- 
tribution function F{t]j\rii) — J^' dip' P(-0'|r/) for all values 
of 77 and ip which constitute the numerical grids for Q{r]) 
and P(V') respectively (discussed in App. |B]). Working 
with the cumulative distribution on a numerical level, as 
opposed to the PDF itself, is far more controlled and cir- 
cumvents divergences in the PDF P{ip), but also in the 
conditional PDF P{il)\rj). The self-consistency condition 
^ is not influenced by this approach. As can be seen by 
inspection, the insulating solution P(e) = (5(e) is always a 
self-consistent solution, equivalent to the ip = Q solution 
in the single site theory. However, in the SF regime, there 
also exists a second, non-trivial self-consistent solution, 
which corresponds to a lower grand canonical potential 
and is therefore the physical solution in this case. Fur- 
thermore, the physical solution is always found to be the 
attractive fixed point of the self-consistency mapping in 
the space of probability distributions, i.e. if the iteration 
procedure is started at any P{'ip) 7^ the successive 

distributions continuously converge towards the physical 
distribution. 

We start the iterative procedure with an initial PDF 
P("'('0), where all the weight is distributed at small, but 
non-zero values of ip, assuring fast convergence in the 
insulating state and in the vicinity of the phase border. 
The distribution in the i-th iteration step for the scaled 
sum of MFPs from the nearest neighboring sites Q^*^ (f?) 
is calculated from P^'-* [ip) using the convolution theorem 
for independent random variables and the FFT algorithm 
(see App. |A|. The new distribution p('+^)(-0) is then 
obtained by integration over 77. Numerically this is done 
by using the trapezoidal rule, as we found that using 
higher order techniques, such as Simpson's rule, is not 
robust and lead to incorrect results if (5-pcaks appear in 
the PDF. 

In the vicinity of the phase border, the computational 
effort increases due to two effects. Firstly, the conver- 
gence is critically slowed down, increasing the required 
number of iterations. Secondly, the discretization of the 
V'-grid plays an ever increasing role and in some param- 
eter regimes directly on the outside of the the MI lobes, 
where the converged distributions are very close to a 5- 
peak at t/i = 0, the numerically determined form of the 
distributions depends on the discretization (resolution), 
which is clearly unphysical. In these cases, we have to 



determine P(V') by examining a sequence of converged 
distributions at ever increasing resolution and define the 
physical distribution as the limit of this sequence. 

III. PHASES OF THE DISORDERED BHM 

A. Phases at T = 

For the disordered BHM, three different phases exist at 
zero temperature: A Mott insulating state, where num- 
ber fluctuations arc suppressed and the particles arc lo- 
calized due to a repulsive interaction. This state exhibits 
a finite energy gap of order U , thus the single-particle 
density of states (DOS) at w = vanishes and the state 
is incompressible. 

If tunncling-induced delocalization dominates, the sys- 
tem is in a condensate (SF) phase, where a macroscopic 
number of particles can lower their energy by condens- 
ing into one single-particle state, thus exhibiting quan- 
tum coherence and leading to a finite condensate frac- 
tion. Within a grand-canonical mean-field description, 
this phase breaks the J7(l)-symmetry of the BH Hamil- 
tonian ([T]) and leads to a non-zero order parameter. The 
phase border from the SF to any of the insulating phases 
is thus determined by SMFT, where finite weight in P{ip) 
moves to finite values of ip. 

Finally, there is the Bose glass phase, where particles 
are localized by an interplay of disorder and interactions. 
However, there exists no single particle state which is 
occupied macroscopically and thus the BG is not a con- 
densate, i.e. the condensate fraction vanishes (/c = 0). 
However, there does exist an extensive number of local- 
ized single particle states, each of which is occupied by 
an arbitrarily large, but not macroscopic number of par- 
ticles. This may be understood as a highly fragmented 
system of incoherent localized 'non-macroscopic quasi- 
condensates'. 

The transition from a BG to a MI is not determined 
from the self-consistent distribution P('0), but by ana- 
lyzing the compressibility n or the single particle DOS. 

We consider the latter quantity within two frameworks: 

1. Considering the purely sing site particle- and 
hole excitations as specified in22. The BG 
extends over the region where P{ip) = 5{ip) 
and values of the chemical potential /i e 
{m -I- e I m g N, e e M,_p(e) > 0}, i.e. the borders 
are independent of J and a direct MI-SF transition 
is possible. 

2. A more detailed analysis presented in2ii22 relies on 
the analysis of an effective Hamiltonian in the sub- 
spaces of localized single particle- and hole excita- 
tions. For finite hopping J in the pure system, these 
hybridize, lifting the degeneracy and form superpo- 
sitions with quantum numbers k. Increasing (de- 
creasing) ^/U , the Mott insulating state remains 
the ground state until the energy difference between 
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the MI state and the fc = particle (hole) state van- 
ishes, at which point particles delocalize and con- 
dense into these states. Now let us return to the dis- 
ordered case: here the local particle (hole) excita- 
tions will not hybridize into fully delocalized states 
with well-defined quasi-monientuni fc, but into in- 
homogeneous states, which depend on the individ- 
ual disorder configuration. It is however possible 
to make exact statements about the eigen-energy 
spectrum for a disordered system in the thermody- 
namic limit: The lowest kinetic energy is obtained 
in locally homogeneous regions and approaches the 
energy of the fc = particle (hole) state of the pure 
system, as the size of this locally homogeneous Lif- 
shitz region increases. Furthermore, scaling pre- 
dicts that the dependence of the kinetic energy on 
the specific boundaries to the Lifshitz region will 
reduce, as the size of the region increases. On 
the other hand, the potential energy of the parti- 
cle (hole) state is minimized, when the local on-site 
energy takes on the lowest (highest) possible value 
over the whole region, i.e. lies at the extrema of 
p(e). Therefore, upper (lower) phase boundary of 
the Mott lobe in the disordered system is obtained 
by shifting the upper (lower) boundary down (up) 
by max({e|p(e) > 0}) (min({e|p(e) > 0})). The MI 
for box disorder in 3 dimension at T = obtained 
by strong coupling theory using this criterion, is 
the area enclosed by the orange phase boundaries 
in Fig. [3l The insulating region outside of these 
lobes, bounded from the SF by the dashed white 
lines, corresponds to the BG. Using this criterion, 
the transition from MI to SF always occurs through 
the BG phase for box disorder within SMFT. 
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FIG. 1: Energy structure of the BHM for large U/J without 
disorder. The levels correspond to the energies of the exact 
many-particle energy eigenstates of the system with interac- 
tions. Due to the [/(l)-symmetry the energy eigenstates can 
be chosen to be of well-defined particle number. For commen- 
surate particle number A^, the ground state (MI) is separated 
from a band of particle- hole excitations by a charge gap, while 
for A'^ -f 1 particles, all eigenstates are delocalized, forming a 
band of hybridized particle states. At T = 0, the single par- 
ticle DOS has finite weight at frequencies corresponding to 



energy differences Ef^^^'^ — E^" ' of transitions between A'^ 
and A'' -I- 1 states. At T > transitions between different ex- 
cited states also contribute, causing the gap to vanish if the 
two bands overlap and the respective matrix elements do not 
vanish. 



p(JV) 



B. Phases at T > 



At finite temperature T > 0, the system is always 
compressible and the incompressible Mott insulator is re- 
placed by a normal (non-superfluid) phase with a ther- 
mally induced compressibility. A central statement in^l 
is, that the disorder-averaged single particle DOS calcu- 
lated by considering purely local excitations only, allows 
for a clear distinction at T > between the BG and 
the normal phase, as it is zero in the latter. However, 
when considering delocalized excitations in the particle 
and hole sector, this statement holds no longer, as there 
are generally degenerate states in the A^-particle particle- 
hole band and the particle (hole) band in the TV -1-1 sector 
(or the A^ — 1 sector for the hole band), which can be seen 
from the single particle DOS in the Lehmann represen- 



Here, E 



(N) 



;th 



denotes the energy of the many- 
^''^ of Hbh, which can always 

be chosen to have well-defined particle number A'^, since 
[HbhtN] = 0, while I is a label for the eigenstate 



particle eigenstate 



in this subspace. 



fi 



(;,fe) 



7/' 



(N+l) 



and 



are the amplitudes for the 

various possible transitions. If the system is in a Mott 
insulating state (considering no disorder for the clari- 
fication of this argument) as shown in Fig. IIII B[ the 
ground state with A'' particles is separated from a band 
of particle- hole excited states by a gap 5{J,U) (corre- 
sponding to the height of the Mott lobe at the respective 
J/U). As transitions between various excited states can 
also occur at finite temperature, the SPDOS thus van- 
ishes if the particle-hole band with N particles, overlaps 
with the particle- or hole band containing A^ H- 1 or A^ — 1 
particles, but the weight is suppressed exponentially by 
S/T. We thus conjecture, that the BG and the normal 
phase are not fundamentally distinguishable, and only 
connected by a crossover. 
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Following along the lines o^, we use the disorder aver- 
aged local single particle density of states (DOS) at zero 
frequency 'p{uj = 0) (see Eq. (jE4p . (jESp ). to determined 
the normal phase / BG crossover at finite T. In the pres- 
ence of disorder the SPDOS containing localized excita- 
tions can be calculated explicitly from the local Green's 
function in this regime (see App. |E|) and leads to 



p(l^,M,A,/3) ^ / dep{e) 



1 



Z(M-e) 

OO 

m=0 

+ mS{io — U{m — 1) + fi — e)] . 

(13) 

Z{^') is the local partition function at an effective chem- 
ical potential n' and the two ^-distributions correspond 
to local particle and hole excitations respectively. Using 
this criterion leads to an overestimation of the MI / nor- 
mal phase over the BG and thereby an absence of the 
BG phase around the tips of the Mott lobes. 

C. Deviations in Finite Size Systems 

Although phase transitions are, strictly speaking, only 
well defined in the thermodynamic limit, crossovers ob- 
served in current experiments may indicate phase bor- 
ders that do not coincide with the borders obtained in 
systems of infinite size. As discussed above and i n^^'^^ , 
the MI/BG phase borders in an infinitely large disordered 
system simply correspond to the shifted phase borders of 
the pure system. This argument relies on the existence 
of arbitrarily large Lifshitz regions, which is clearly no 
longer given for finite systems. The phase diagram for a 
finite system therefore strongly depends on the specific 
disorder realization, with the critical values for the phase 
borders becoming random variables. Therefore, a better 
question to ask for a finite system for instance is: For a 
randomly chosen disorder realization in a system consist- 
ing of L sites, what is the probability Pg that the energy 
gap will be lower than a given value D? In the limit of 
J/U ^ this probability is given by 



PgigaiX D) = I - {1 - Pd)"- 
with the restriction that < D < U/2 and 

OO pmU+D 

Pd=^ d^'p(/i- /i'). 



(14) 



(15) 



With increasing J/U the MI/BG phase border in the 
J/[/-/i/C/-diagram remains a random variable, but with 
a reduced steepness in the slope for finite systems, as 
shown in2ii2^. The BG/SF would also be very likely (in 
a statistical sense) to move to larger critical values of J 
in a finite system, as the occurrence probability of 'rare 
events', favoring a SF, is suppressed. 



IV. BOX DISORDER 

Results for box disorder at T = obtained by SMFT 
have been presented in^^, here we extend the phase di- 
agram by taking collective excitation in the BG into 
account and focus on finite temperature. All numer- 
ical results presented in this paper are for a three- 
dimensional cubic lattice (Z=6). Box disorder is charac- 
terized by a constant probability density for the on-site 
energies p{e) = ^6(A/2 — |e|) over a bounded inter- 
val of width A, where Q{x) denotes the Heaviside func- 
tion. Changing the disorder strength may have an infiu- 
ence on the system properties on various levels. Com- 
ing from the insulating state within a local mean-field 
picture, where the total state of the system is a direct 
product of local Fock states minimizing the total en- 
ergy, an increase in disorder reduces the smallest pos- 
sible particle or hole excitation energy. The local en- 
ergy gap for a site with an effective chemical potential fi' 
and a Fock state as ground state \n = max([/i'/C/] , 0)) 
is iSparticic = max([/z'//7] , 0) [/ — fi' for adding a parti- 
cle and i?i„i„ ~ fi' — {ma.x{\fi' /U~\,0) — 1)U for creating 
a hole, if /i' > 0. In the presence of disorder in an in- 
finitely large system, this implies that the energy gap is 
necessarily reduced by A/ 2 and vanishes in the J — >■ 
limit as soon as the interval of realizable effective chem- 
ical potentials contains a positive integer multiple of U, 
i.e. N n [(^i - A/2)/C/, {fi + A/2)/U] 7^ 0. 

Within the picture of purely local excitations, this ar- 
gument is independent of temperature and directly re- 



KU for T=o.ooi U 
KU for T=o.02 U 
KU for T=o.2 U 
DOS p(to=o) for T=0.2 U 




A (units of U) 



FIG. 2: (Color online) Compressibility k and local single par- 
ticle DOS in an insulator for box disorder at ^ = I.IU at 
various temperatures (in units of U^^). At zero tempera- 
ture K vanishes in the MI, but not in the BG and may thus 
be taken as a quantity to distinguish between them. At fi- 
nite temperature however, k becomes non-zero in the normal 
phase and the sharp features become rounded off by thermal 
fluctuations. In this case, the disorder averaged DOS of lo- 
cal excitations p(i^ = 0), which retains all sharp features at 
T > 0, can be used to determine the crossover from a BG to 
the normal phase, where it remains zero. 
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fleeted by the disorder averaged local DOS 



p(^,M,A,/3) = -Y^{m + l)^ ^ 

m=0 



Z[Um — ui) 



(16) 



In the case of bounded disorder, this directly implies 
that the system cannot be in a gapped state as soon as 
the carrier width (A for box disorder) equals or exceeds 
U and that the system is then always in the BG or SF 
state, i.e. the MI and normal state can only exist at 
A < U (and not at all for unbounded disorder). When 
the system is either of these two insulating states at zero 
temperature, the compressibility K = §2 is determined 
by ^(0; = 0) within a picture of local excitations, where 
locally degenerate sites change their occupation number 
by ±1 if yu is slightly altered. Since for finite temperature 
thermal particle or hole excitations are also present when 
the system is gapped, k is always driven to a finite value. 
The behavior of k and p(a; = 0) is shown in Fig. [2] at 
the specific chemical potential of /i = I.IU. Starting at 
A = 0, the system is in a gapped MI / normal state, but 
increasing A, the system undergoes a phase transition at 
a critical disorder strength of A = 0.2[/. 

At this point the local DOS p(w = 0) takes on a non- 
zero value, as two different local Fock states can become 
degenerate at the edge of the box disorder distribution. 
At low temperatures the compressibility exhibits a sig- 
nificant change from an exponentially small value in T in 
the normal phase to a large value due to the existence of 
local configurations of on-site energies which give rise to 
almost degenerate states with different particle number. 
There, a small change in the chemical potential leads to 
a local jump in particle number. In the zero temperature 
limit the compressibility locally behaves as k oc ^, since 
with increasing disorder strength the statistical weight 
of these events is reduced. When the effective chemi- 
cal potentials at the borders of the probability distribu- 
tion enter a region corresponding to a new local parti- 
cle number ± A^ the compressibility increases with a 
jump, corresponding to subsidiary phase transitions be- 
tween different BG phases at T = 0, which turn into less 
pronounced crossovers with increasing T. Considering 
the limit of strong disorder A/U — 00 at fixed ^/U in 
Fig. [21 the system approaches a state in which only half 
the number of sites are occupied and limA^oo — 

However, disorder does not only affect local excitation 
properties, but leads to an intricate interplay between 
the hopping and interaction energy scale, influencing the 
overall coherence properties of the system. Whereas an 
increase in temperature or interaction energy generally 
tends to counteract the formation of a condensate, this 
is true in most, but not all scenarios when increasing the 
disorder strength. In certain parameter regimes at suf- 
ficiently low temperature, an increase in A can actually 
lead to the formation and stabilization of a condensate 
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FIG. 3: (Color online) Finite temperature phase diagrams for 
box disorder with A — 0.5U for three temperature regimes: 
Zero temperature (upper row), intermediate T = 0.03U (mid- 
dle row) and high temperature regime T = 0.2U (lower 
row). The left column shows the expectation value of the 
self-consistently determined distribution, which is the order 
parameter for the SF-insulator transition. It vanishes in the 
MI, BG, and the T > normal phase. In the T = diagrams, 
the MI/BG phase borders from strong coupling theory are in- 
dicated by the orange solid lines. The BG phase is marked by 
white stripes, at finite T > this is determined by a vanish- 
ing local DOS only (i.e. not by strong coupling theory). At 
finite temperature there is only a crossover between the nor- 
mal and BG phase, and the borders determined by the DOS 
of purely local excitations is marked by the horizontal white 
lines. The right column shows the compressibiUty k in units 
of C/~^. Whereas this is a suitable quantity to distinguish the 
MI (incompressible) from the BG at T = 0, it is non-zero, but 
exponentially small in T in the MI in the intermediate regime 
and of order U^^ in the high temperature regime. The white 
borders indicate the transition to the SF, determined from 
the respective diagram on the left. 



(i.e. disorder induced condensation), as previously pre- 
dicted by various methods, including SMFT— . 

In Fig. [21 the effect of increasing temperature is exem- 
plified in three phase diagrams and the compressibility 
in the three main regimes: the zero, low and high tem- 
perature regime at fixed disorder strength. The tips of 
the MI /normal lobes remain almost unchanged under 
an increase in T, while the BG region between the initial 
lobes is strongly enhanced and stabilized, even at a very 
low temperature of kT = 0.03U (central plots of Fig. [3]). 
Since the SF/insulator phase border still possesses a dis- 
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tinctive lobe structure, but the BG region is strongly 
enhanced in contrast to the T=0 case (upper plots of 
Fig. [3]), wc therefore refer to this situation as the inter- 
mediate temperature regime. Furthermore, a large value 
of compressibility is still a strong indicator for the BG 
at this temperature, suggesting that it may still useful as 
an indicator of the transition between the BG and the 
normal phase in experiment. In the high temperature 
regime kT ^ 0.2U, on the other hand, (see lower dia- 
grams in Fig. |3]) the compressibility is large throughout, 
approaching the value k — >■ in the high T limit. The 
typical temperature of T « 0.2U at which the MI and 
normal phase melts is consistent with previous previous 
studies^. In this temperature regime the lobe structure 
of the SF/insulator phase border is totally wiped out and 
the critical border to the SF follows a {JZ)c (x fi~^ decay 
with increasing filling. 

The SF/BG transition is highly sensitive to the system 
size. When the system approaches the SF phase from the 
BG phase, the localized single particle orbitals occupied 
by a large number of bosons increase in size and are oc- 
cupied by an ever increasing, but never extensive number 
of particles in the BG phase. At the transition point the 
localization length, being a measure for the size of these 
orbitals, diverges and driven by percolation, phase coher- 
ence between neighboring orbitals is established, eventu- 
ally driving the system into the SF phase. In finite size 
systems, the detection of the transition point thus criti- 
cally depends on the system size, as the BG phase may 
be mistaken for the SF phase if the localization length 
is larger than the system size. SMFT has the advantage 
that it is constructed in the thermodynamic limit for an 
infinitely large system in the grand canonical ensemble 
and takes all possible disorder realizations into account 
within a functional description for the probability distri- 
butions. The numerical error in the discretization per- 
formed for the distribution P{ip) is well controlled and 
to be distinguished from the finite size deviations made 
in real space calculations. 

In the phase diagram for box disorder Fig. [31 SMFT 
does not give rise to a direct transition from the MI to 
the SF, if the extended criterion for the MI/BG border, 
including collective excitations in rare regions is used. 
This is furthermore demonstrated in the finite tempera- 
ture phase diagram at constant filling n = 1 in Fig. 2] 
Due to the absence of a clear distinction between the 
normal and BG phase at T > 0, the orange border only 
indicates a crossover between these regimes, but would 
go over into a MI/BG phase border for T = 0. At any 
A > a finite BG region intervenes between the MI and 
SF phases. 

The question whether this transition always occurs via 
the BG phase has been a highly debated topic since the 
introduction of the disordered BHMi, and was estab- 
lished for the one and two dimensional case^^ii^i^. In 
a recent work^i it was shown that this scenario is true 
in any finite dimension for bounded disorder, due to the 
statistical certainty that any possible configuration of on- 



site energies for a cluster of sites will occur in the limit of 
an infinitely large system. In the previous work^ these 
collective excitations in the BG phase were not consid- 
ered and within a simpler framework of purely local parti- 
cle and hole excitations, a direct transition was predicted. 

It is known that arithmetically averaged MFT, as well 
as SMFT providing an improvement in any finite dimen- 
sion, both become exact in the limit of infinite dimen- 
sions, where no BG exists at T = 0. However, as argued 
recently^ the theorem of inclusions guarantees the ex- 
istence of an intervening BG phase between the SF and 
the MI in any arbitrarily high, but finite number of di- 
mensions for bounded disorder. It is instructive to un- 
derstand the decrease of the BG region, including the 
collective excitation in the Lifshitz regions in terms of 
percolation physics: For any finite dimension the outer 
border between the BG and SF phase specifies the criti- 
cal value ( J/?7)„it.i specifies the lowest energy at which it 
becomes energetically favorable for the particles to form 
a global condensate (in SMFT this is the border where ip 
takes on a finite value). The border between the MI and 
BG inside the global insulator, specifies the critical value 
(J/C/)„it.2 at which it becomes possible for the system 
to form large local superfluid patches (locally resembling 
pure systems) without phase coherence between differ- 
ent patches. With increasing dimensionality of the sys- 
tem, the connectivity between different patches increases 
(percolation is enhanced) and the required tunneling en- 
ergy ( J/[/)crit.i to form one large percolated patch, i.e. a 
global condensate, decreases. In the limit of high di- 
mensions, the critical values of J at which these two 
phenomena occur approach the arithmetically averaged 
mean-field valued 



JZci^l) =A 



n In 



+ {n+1) In 



1 - ?i + /i + A/2 
1 - ?i + ^ - A/2 
- H + A/2 



n- fi- A/2 



(17) 



and the BG disappears, where n is the filling and fi and 
A are given in units of U. 

In Fig.|3]we present a phase diagram calculated at fixed 
density n = 1 in the low temperature regime T — 0.03U. 
At every point in the diagram, the self-consistent dis- 
tribution is calculated for a fixed /i, enabling the calcu- 
lation of the density {n{fi, A,U, J,T)) . Thereafter /i is 
iteratively determined using Bidder's algorithm^ until 
the density obtained from SMFT does not deviate more 
than An = 0.005 from the specified density. In Fig. 01 the 
disorder averaged MFP = J tp P{^)dip (within SMFT, 
this is exactly v7^) clearly shows the usual SF/insulator 
phase transition (at fixed low temperature) along the line 
A = 0, where the disorder localizes the particles with in- 
creasing U/J. Moving outwards into the A/ J at fixed 
interaction U/J, the condensed phase is surprisingly ro- 
bust, surviving local on-site fiuctuations A several hun- 
dred times larger than the hopping energy J, as pointed 
out in a recent work2^. This can be understood from 
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the bosons filling up the low-lying sites and forming a 
'background sea' via the repulsive interactions, creating 
an effective smoother potential in which it is energeti- 
cally favorable for the remaining bosons to delocalize^. 
A remarkable effect at sufficiently low temperature is the 
appearance of a SF lobe, protruding into the insulating 
domain at finite A. In this regime the interplay between 
disorder and interactions is non-monotonic in these two 
effects, and, for the regime AO ^U/ J ^ 85, an increase in 
A drives the system into the SF phase, delocalizing the 
particles. This effect can be understood from the pure 
BHM n/U-J/U phase diagram and relies on the existence 
of a lobe structure, i.e. requires a sufficiently low T. To 
keep the particle number constant with increasing A, fi 
is required to increase. In certain regimes the majority of 
sites in the system may enter from an insulating regime 
between the lobes into a regime inside the lobes, thereby 
favoring condensation. Qualitatively, the SMFT phase 
diagram agrees well and shows the disorder induced SF 
lobe, as found in recent QMC calculation^^ at T on 
relatively small lattices (L = 8 x 8 x 8). 

At large A the order parameter is non-monotonic in 
U/ J, vanishing at sufficiently small U/ J which indicates 
a transition into an Anderson localized state, where the 
localization almost exclusively disorder-induced. How- 
ever, the region of extremely small U/ J is problematic 




FIG. 4: (Color online) Phase diagram for box disorder at fixed 
density n = 1 in the U/J-A/J plane at T = 0.03(7 showing 
the mean order parameter if; = \pfc- Reentrant superfluidity 
is reflected by the protruding SF lobe (we also find a second 
less pronounced lobe at higher A), where increasing A can 
drive the system through a sequence of SF-insulator transi- 
tions. The orange line specifies the crossover from the normal 
to the BG phase (at T = this line becomes the MI-BG phase 
border) as determined by by shifted mean-field phase borders 
at this temperature, giving a better approximation than the 
strong-coupling approach for three dimensions. The simpler 
criterion of looking at purely local particle and hole excita- 
tions Ea.ll6lwould lead to the MI/BG phase border at f/ = A 
and predict a direct MI-SF transition at small A. Including 
collective excitations in the BG, as done here, always leads 
to an intermediate BG phase between the MI and SF in the 
T = Umit. 



as A/ [7 and [ijU diverge, since very few sites have to 
contain an ever increasing number of particles to keep 
the disorder-averaged density fixed, when asymptotically 
half the number sites (due to the symmetry of the box 
distribution p(e)) have such a high effective on-site en- 
ergy, that they contain no particle. Due to the diverg- 
ing local occupation number, this limit transcends the 
constraints imposed in the derivation of the BHM in an 
optical lattice and is, in this sense, unphysical. 

V. SPECKLE DISORDER 

Although a homogeneous box distribution is most com- 
monly used for disorder calculations in theory, it is cur- 
rently not an experimentally feasible choice. In this sec- 
tion we discuss and compare the results for a realistic dis- 
order distribution created by a speckle laser to those of a 
box disorder distribution. A laser passing through an in- 
homogeneous disordered plate leads to a disordered opti- 
cal potential, which is the Fourier transform of the disor- 
dered pattern on the plate. In recent experiments, it has 
become possible to reduce the autocorrelation length of 
this disordered potential to the order of the lattice spac- 
ing (< With this experimental achievement, the 
priorly most criticized artifact of a speckle laser for cre- 
ating uncorrelated disorder has been overcome, thereby 
making speckle potentials the most promising method 
for future disorder experiments in optical lattices. The 
resulting distribution for uncorrelated on-site energies is 
well approximated by 

p(e) = ^e-/^ (18) 

Although it may be argued that an optical speckle po- 
tential in experiment is fundamentally bounded by its fi- 
nite size, it is only essential that the width of the on-site 
energy distribution exceeds U (which is fulfilled in essen- 
tially all experimentally relevant regimes). This justifies 
the use of (|T5)) . 

To treat this disorder distribution using SMFT, it 
is useful to perform a transformation of variables 
x(e) = — e^^/'^, which on a formal level transforms the 
SMFT conditional probability functions into a form anal- 
ogous to homogeneous disorder. This step enters only on 
the level of calculating the conditional cumulative distri- 
bution function (CDF), 

F(ip\r]) = \im[ dx Q(ip ~ gU + A \n(-x),ri)). (19) 

Apart from this, the SMFT method remains identical to 
the homogeneous disorder case. Similarly, arbitrary dis- 
order distributions may also be incorporated into SMFT, 
although an analytical transformation of the random 
variable will not exist in general. 

In contrast to box disorder, which has been the dis- 
tribution primarily focused on so far when considering 
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FIG. 5: (Color online) Typical self-consistent distributions 
P('i/') for speckle disorder for A = If/, = 1.2C/, Z = 6, and 
= 0.3C/ at different temperatures. For low temperatures 
the mean-field parameters in the condensed phase are robust 
against finite temperature fluctuations, but with increasing 
temperature the system is eventually driven into an insulating 
BG phase, as indicated by the distribution at P(V') = <5(7/;) in 
the at T = 0.35[/. 
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the disordered BHM, speckle disorder is unbounded and 
arbitrarily high values of e have a finite probability to oc- 
cur. This leads to the effect (strictly only possible for an 
infinitely large system) that turning on the disorder by 
an arbitrarily small amount immediately changes a large 
part of the phase diagram from the MI / normal phase 
into the BG phase. 

In a physical picture, the effective chemical potential [il 
can then take on integer multiple values of \J (the on-site 
potential can become arbitrarily high) with a non-zero 
probability density, where local Fock states with different 
particle number become degenerate, leading to a finite 
compressibility and local DOS /3(w = 0) at T = 0. 

In the absence of disorder and at sufficiently high tun- 
neling coupling J, a macroscopic number of particles oc- 
cupy the |fc = 0) Bloch state. Within an effective, sym- 
metry breaking Gutzwillcr description the local order pa- 
rameters 1^1 = (6); then take on a finite and constant 
value, reflecting the translational symmetry of the sys- 
tem. Within SMFT this state is characterized by a 8- 
distribution P('!/') = ^('!/' ~ V'o)) where '00 is the order pa- 
rameter of conventional bosonic Gutzwiller theory. Turn- 
ing on the disorder in such a system in the SF state breaks 
the translational symmetry of the system, i.e. the con- 
densate state deviates from the fc = Bloch state, which 
is reflected by the distribution of MFPs P{i)) taking on 
a flnite width. Initially for weak disorder, an increase 
in disorder always leads to a broadening of P{ip), but 
for stronger values of A the system may eventually be 
driven toward an insulating state, driving P{ip) — >■ S^ip) 
and thereby decreasing the fluctuations in the MFPs. On 
the other hand, increasing the temperature suppresses 
the SF and leads to a decrease of the MFPs above a cer- 



FIG. 6: (Color online) Disorder averaged MFP t/) = 
J d'ip ^j) P{ip) characterizing the SF-insulator transition (left 
column), and compressibility (right column, in units of U~^) 
in the JZ/U-f-i/U-plane for speckle disorder. Diagrams are 
shown for increasing disorder strength (A = 0.1(7, A = 0.3U, 
A = If/) in the zero or low temperature regime (T = 0, 
T = 0.05f/, T = 0.05f/) for the (upper, middle, lower) row 
respectively. In contrast to the behavior for box disorder, 
the structure of each lobe does not change symmetrically, but 
rapidly extends in the direction of decreasing /i/t/ with in- 
creasing disorder strength A. Since the MI / normal phases 
do not exist for speckle disorder, the insulating region (black 
in the left figures) is always a BG. The white lines denote the 
SF / insulator phase boundaries, indicating where ip takes on 
a finite value. 



tain temperature, up to which the SF remains stable, as 
shown in Fig. [S] 

The influence which speckle disorder has on the n/U- 
J/ [/-phase diagrams is shown in Fig. |6l In contrast to 
box disorder, where the distribution of on-site energies is 
symmetric around n and the insulating lobes give way to 
the SF in the same way on the upper and the lower side 
of the lobe, the insulator forms on the lower side of the 
lobes with increasing A for speckle disorder. This can be 
understood from the fact that only lower values of the 
effective local chemical potential can occur. 

For strong disorder A ^ [/, the lobe structure of the 
insulator / SF phase boundary is washed out, which is 
similar to the effect of finite temperature. For speckle 
disorder, n cannot be used to identify a phase transition, 
since it is non-zero in both the BG and the SF. A ques- 
tion of interest, regarding the disappearance of the MI 
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/normal phase for an arbitrarily small amount of speckle 
disorder, is how the compressibility behaves as a function 
of A, as some interpretation is needed, that an 'infinites- 
imal amount of disorder' can instantaneously convert the 
whole Ml/normal area of the phase diagram into a BG. 

In Fig. [71 k(A) is shown for different parameters to 
clarify this dependence. In the insulating state, the com- 
pressibility can be calculated explicitly (|D3|) and one ob- 
tains 
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Essentially, two different scenarios have to be considered. 
First, if /i/C/ is positive and integer, k diverges in the 
limit of vanishing disorder, as is well known from the 
pure BHM phase diagram, where the density is a step 
function in yujlJ for J = 0. This is equivalent to the 
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FIG. 7: (Color online) The compressibility, local single par- 
ticle DOS and chemical potential as a function of disorder 
strength A for speckle disorder at constant filling. Compari- 
son of the upper two plots shows a diverging compressibility 
and local DOS for A — >■ for non-integer filling (the sys- 
tem remains superfiuid), while for integer filling it drops to 
an exponentially small value in A at a value of A ~ 0.1J7 
and vanishes in this limit. The lowest figure shows the same 
quantities at finite temperature kT = O.Sf/, where thermal 
fluctuations have totally smeared out the sharp features in 
the compressibility. However, these persist in the local DOS, 
although their position is changed due to a temperature in- 
duced shift in ^. 



case for non-integer, fixed particle density n. where the 
system remains SF for any non-zero J. Second, if [ijU is 
non-integer, the compressibility vanishes with decreasing 
A, as the system approaches a point in a Mott lobe away 
from the border. This corresponds to the case of fixed, 
integer-valued density n. 

We will now discuss the behavior of k and p(a; = 0) 
at fixed particle density n, shown in Fig. [71 Keeping 
the density constant with rising disorder, requires the 
chemical potential to be increased, as an ever increasing 
number of sites shifts to weights with lower occupation 
numbers. At every point when [ijU passes a positive 
integer number, a new Fock state becomes potentially 
occupied, but with an ever decreasing statistical weight 
as A increases. As a result the compressibility experi- 
ences a jump at each of these points, as highlighted by 
the gray regions in Fig. [71 where lijU (dotted gray lines 
in Fig. [7]) passes an integer value. This leads to the char- 
acteristic series of ever smaller kinks in k(A). At finite 
temperature, these features in n are smeared out over a 
typical scale of fcT, whereas the sharp features in the lo- 
cal single particle DOS survive at T > within SMFT 
(solid blue lines in Fig. [7]). 




FIG. 8: The disorder averaged single particle density of states 
p(ijj = 0, \i) at zero frequency for T = 0.05(7 for different 
speckle disorder intensities. For A = this consists of a 
sequence of (5-peaks at positive integer values of iijU , however 
for any A > this quantity is non-zero for > and the 
system is in the BG phase. 

To clarify the effect speckle disorder has on p(cj = 0) 
and the immediate disappearance of the MI / normal 
phase at any A > 0, the local DOS is plotted for 
weak (A — O.IC/), intermediate (A = 0.5(7) and strong 
(A = 2(7) in Fig. [H as a function of /i. In the pure sys- 
tem p(a; = 0, [I) consists of a sum of (5-peaks at integer 
values of /i/(7, i.e. at these values of the chemical poten- 
tial there are two degenerate Fock states |n = [xjU) and 
\n = [ijU -1-1) at all sites in the insulator and the local 
single particle DOS diverges. As soon as speckle disorder 
is turned on, these (5-peaks are broadened according to 
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the on-site energy distribution (|T8|) and p(cj = 0) takes 
on the form of a sequence of superimposed exponential 
functions, each decaying with the constant A. Prom this 
it is clear that p(a; = 0) takes on a non-zero value as soon 
as A > at any ^ > 0, although it is exponentially sup- 
pressed for most values of ^ at weak disorder A ^U. At 
zero temperature the different amplitudes of the various 
peaks in 'p{uj — 0, fi) at integer fi can be exclusively at- 
tributed to the ^/n factor from the action of the bosonic 
operators, whereas at T > the amplitudes (but not the 
positions of the sharp features) may also be modified by 
the Boltzmann factors in (IE4I). 
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FIG. 9: (Color online) Main image: Phase diagram in the 
T/ J-C// J-plane for fixed A — 10 J at fixed filling n = 1 in 
three dimensions. For this specific disorder strength the SF 
region is enlarged by the disorder in comparison to the pure 
case (i.e. disorder induced condensation, see Sec. IV A")) . In the 
lower left corner, the SF exists at sufiiciently low T and A and 
the value of disorder-averaged tp is color-coded. Outside the 
SF region, the system is always in a BG phase, but undergoes 
two crossovers from a regime with exponentially low com- 
pressibility kU for intermediately low 0.02 ^ T/U ^ 0.065 
(for A/J = 10) into a strongly compressible BG regime, in 
the limits of both high and very low T/U. To clarify the 
quantitative behavior, the compressibility k, the chemical po- 
tential fj. and the local single particle DOS p{uj — 0) along the 
dashed line U/J = 100 is shown in the inset. 

Fig. ini shows the phase diagram at constant disorder 
strength A/J~l and constant density n = 1 in units 
of J. The SF region prevails in the lower left region 
at low temperature and weak interaction U/J. In this 
parameter regime the disorder stabilizes the SF phase, 
actually extending the SF region of the phase diagram 
in contrast to the pure (A = 0) case (disorder-induced 
condensation, see Sec. IV Ap . All of the non-SF region 
in Fig. [S] is a BG, since we are dealing with unbounded 
disorder with A > 0, but we can identify a weakly and 



two strongly compressible regimes in the phase diagram. 
Since P{tp) = S^tp) in the BG, the energy scale J can- 
not influence thermodynamic quantities beyond a scaling 
relation, implying that the compressibility kU may be a 
function oi T/U only. Therefore the compressibility has 
a radial structure and it suffices to consider its behavior 
along a single line (such as U/J= 100, depicted in the 
inset of Fig. in]). This reveals that there are three regimes 
in this phase diagram: at high temperature T/U ^ 0.065 
the system is strongly compressible {kU is of order unity), 
as thermal fluctuations have wiped out the sharp peaks 
over a wide range in ^. At intermediately low tempera- 
ture, the compressibility is exponentially low, as the typ- 
ical thermal excitation energy scale does not suffice to 
excite the majority of sites to higher states. Somewhat 
surprisingly, within a certain parameter regime for A and 
the density n, the system undergoes a further crossover 
at very low temperature T/U ~ 0.02 into a second highly 
compressible regime. At T = and at integer filling, the 
chemical potential fi/U approaches an integer value from 
below (green dotted line in the inset of Fig.[S]). At these 
points K diverges in the pure limit of A 0, which re- 
veals that the compressibility is grows with the inverse 
of A and only persists in the limit T — > if the density 
n is integer. 



A. Disorder-induced reentrant Superfluidity 

At low temperature and fixed density, we find that 
an increase in A can actually drive the system from an 
insulating into a condensed state within a certain window 
of U /J, as shown in Fig. [TO] and Fig. [TT] Within a very 
small window oiU / J , the system may be driven through 
an additional sequence of BG and SF phases. 
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FIG. 10: (Color online) Low temperature kT = 0.03(7 phase 
diagram showing the disorder-averaged MFP -0 for speckle 
disorder at constant filling n = 1 in three dimensions. Multi- 
ple reentrant behavior can be seen within a small window of 
U/J. The red line at A = indicates the presence of MI / 
normal state, for any A > the insulator is a BG. 

On the A = line, the usual SF - insulator transition 
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occurs, whereas along A for small values of (7, disorder 
only suppresses the SF slightly up to reasonably large val- 
ues of A/J. For even larger disorder values, the system 
will eventually undergo the transition into an 




U/J 

FIG. 11: (Color online) Compressibility k in units of 1/U at 
kT = 0.03C/ for speckle disorder at constant filling n = 1 m 
three dimensions (same parameters as in Fig. llOp . The in- 
compressible normal phase only exists on the line A/J = 0, 
but the region where k is exponentially small in the insulator 
extends linearly with U/J. The white line is the SF - insula- 
tor phase border, where %p becomes zero. The compressibility 
along the blue dashed line can be understood as the low tem- 
perature case and compared to k in the lower two plots of 
Fig. [T] where the sharp features have been washed out by 
temperature, but are still pronounced peaks. 

Anderson localized state. In comparison to the corre- 
sponding phase diagram for box disorder Fig. 01 the ex- 
truding superfluid lobe appears at a considerably smaller 
values of A w 40J for speckle disorder (A « lOOJ for 
box disorder, not the different A-axis scale). However, it 
should be noted that the measure A is not the same (in a 
statistical sense) for the two disorder types: whereas for 
speckle disorder A is also the standard deviation (std), 
in the case of box disorder the std is only « 0.29 A, 
corresponding to weaker disorder for the same A. Ex- 
amining the position of the superfluid lobe in units of 
the standard deviation, actually reveals that it appears 
at slightly smaller disorder strengths for box disorder. 

In the vicinity of the phase border, k is generally larger 
in the SF state than in the neighboring insulator, as a 
shift in /i leads to a locally continuous shift in density. On 
the lower side of the lobe at disorder strengths A ^ 40 J, 
the compressibility jumps across the phase border (i.e. 
second order transition). At higher values of A, k is al- 
most unaffected by the disappearance of the condensate, 
i.e. the compressibility is almost exclusively induced by 
disorder. In the insulating state the ensemble of mean- 
field states determined by SMFT cannot depend on J 
(since different sites are only coupled via ip). Therefore 
any physical quantity, such as k, can only be a function 




FIG. 12: (Color online) Main background image: same ip 
phase diagram as in Fig. llOl but in the experimentally relevant 
high temperature regime kT = 0.3(7 at n = 1, where disorder 
cannot induce condensation. Insets: disorder averaged MFP 
(a), compressibility (b) and local DOS p(tJ = 0) in the ^i/U- 
JZ/U-plane at the same temperature T — 0.3U for A = 117. 
In this regime most interesting structure has been washed out 
by thermal fluctuations and disorder and an increase in either 
f/, A or T always suppresses condensation. 

oiT/U in this regime (reflected in the radial structure of 
K in the insulator). 

In the high temperature regime relevant for current ex- 
periments (see App. lVl"51 for an approximation of T), the 
lobe in the A/J — C//J— plane at fixed density has van- 
ished completely and an increase in disorder always coun- 
teracts the condensate formation. This might explain the 
recent finding, that no reentrant behavior or disorder- 
induced condensation was observed in experiments so 
far—. Furthermore, we calculated the /i/C/ — J/f/— phase 
diagram for this temperature and disorder regime, where 
the lobe structure is also fully washed out and the sys- 
tem is dominated by thermal fluctuations, as seen in the 
insets of Fig. [T^ We therefore conclude that an upper 
critical temperature, which may depend on the filling, 
exists for the occurrence of reentrant superfiuidity. Our 
temperature estimation suggests that T is too high in 
current experiments to observe this effect, which agrees 
with recent experimental evidence^. 



B. Hopping Disorder 

In addition to diagonal on-site energy disorder, it is 
possible to incorporate off-diagonal hopping disorder— 
into the SMFT formalism. In this case, the local hop- 
ping energy J{ij) between site i and j becomes a random 
variable, each described by a distribution pj{J), which 
we assume to be independent of the on-site energies e^. 
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This leads to the BH Hamiltonian with on-site, as well 
as hopping disorder 

^BH = -^J(.,,)(6^6j-+h.c.) 

+ J2ie^-^i)n, + -Y,n^in^~l)■ (21) 

i i 

The PDF pj{J) has been calculated for a speckle disorder 
potentia l^^'^^ using imaginary time evolution. 

In contrast to diagonal disorder (on-site energy or in- 
teraction), where the fluctuations are incorporated into 
the conditional PDF (jSj before the iteration procedure 
of the SMFT self-consistency equations, the hopping dis- 
order acts as an additional source of fluctuations dur- 
ing this iteration procedure, methodically entering at a 
different point in the method. The corresponding mean- 
field Hamihonian "H*""' , (analogous to ©) now depends 
on the rest of the system only through the new random 
variable 

z 

where both Jj and tpj are random variables, the latter 
again being assumed to be distributed according to P{ip)- 
It is therefore convenient to introduce the an intermedi- 
ate random variable (j) = J^^j distributed according to 
the PDF P^{(j)). As explained in App. O this new PDF 
can be expressed explicitly in terms of the PDF's P{ip) 
and p. J (J) as 

PM= j dxpj{enP{^-e-n, (23) 

numerically allowing the use of the FFT algorithm on a 
suitable grid. Once P^iiip) is known, the PDF for the ran- 
dom variable 77 = X^z^i "^^^ subsequently be calculated 

by 

QW = ^ jdtmf e-^'\ (24) 

where 

d{t) = j d<l,P^{cp)e:'t (25) 

is the characteristic function of P^ [(j)) . For the numerical 
computation of the previous two, the FFT algorithm can 
be used. 

The conditional PDF ([5]) P{ijj\rj)^ incorporating the ef- 
fect of any diagonal (here: on-site) disorder remains un- 
changed under the inclusion of hopping disorder and is 
calculated before the SMFT iteration procedure in full 
analogy to the previous case for the relevant on-site dis- 
order type p(e). The final self-consistency equation, clos- 
ing the iteration procedure is also left unchanged to the 
previous case © 

P{i,) = f dr/Q(7?)P(^|r;), (26) 



except that the Q{r]) entering is calculated from (|24)) . 

Of course the case of pure on-site disorder with a con- 
stant hopping energy Jq can be obtained as a limit of this 
extension by setting pj{J) = d{J ~ Jo), causing and 
(Pi)) to reduce to 

The numerical iteration procedure is carried out anal- 
ogously to Sec. Ill A[ except that only a single equidistant 
grid, as restricted by (P^ . can be used. This limits the 
numerically obtained precision of P{ip) at small values of 

In our calculation we use the appropriately scaled dis- 
tribution pj{J), matching the on-site disorder strength, 
as obtained by Zhou and Ceperley^ for a lattice depth 
of s = 14i?/j and at a speckle disorder strength of 
So = lEn, where we assume that the standard deviation 
in pj{J) is proportional to the disorder intensity sd (as 
motivated by their analysis). Since a change in the dis- 
order strength does not change the most probable value 
of the distribution^ Pj{J), wc model the hopping disor- 
der distribution to have its most likely value at the given 
value J/U, with the width being independent thereof. If 
any weight lies at values of negative J, this is set to zero 
and the distribution is subsequently renormalized, how- 
ever this is only relevant for points deep in the insulator 
and does therefore not influence the results in any way. 
However, this only occurs at very small values of J/U 
deep in the insulator, making this formal alteration of 
Pj{J) irrelevant to the result. 

As shown in the low temperature phase diagram in 
the J/C/-^/J7-plane fFigfl^. the additional inclusion of 
hopping disorder leads to a stabilization of the SF phase 
and a small shift in the phase boundary. For disorder 
distributions in the typical experimentally relevant pa- 
rameter regime, the standard deviation of the hopping 
parameter distribution pj{J) is three orders of magni- 
tude smaller than the standard deviation of the on-site 
energy distribution p(e) (as found in^, specifically for 
the distributions used here (Jj/cr^ = 0.0014). This ex- 
plains the minor, but clearly resolved modification of the 
phase boundaries in contrast to the pure on-site disor- 
dered case, shown in the upper right inset of Fig. [T31 

VI. INCORPORATING EXPERIMENTAL 
ASPECTS 

A. Temperature estimation in an optical lattice 

The initial temperature in the trap, prior to the optical 
lattice ramp-up, can be determined from the expansion 
profile of the cloud. If the ramp-up of the optical lattice 
is performed adiabatically, the entropy of the system is 
conserved and can actually lead to a cooling of the atoms. 
The initial entropy of a weakly interacting cloud in the 
trap using Bogoliubov theory leads to the expression^ 

= kB E (J^^ HI e^^p]) . (27) 
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B. LDA incorporating trap effects 

To compare the results obtained via SMFT to exper- 
imental data on a quantitative level, we performed a 
LDA+SMFT calculation to incorporate the effect of the 
trapping potential. Two effects of the lattice laser beams 
arc taken into account: The red-shifted lattice laser beam 
with a Gaussian profile and beam width wq leads to an 
attractive potential via the ac Stark effect. Furthermore, 
the local energies of states within a localized Wannier 
basis are also increased in regions of high intensity, lead- 
ing to the renormalized effective lower trapping frequency 
within an harmonic approximation^. 



§(2.-Vi). 



-1 o 



(30) 



With the addition of an external magnetic trap with 
trapping frequency w^j^g, the total trapping frequency, 
which the atoms are exposed to is given by 



FIG. 13: (Color online) Diagrams showing the effect of hop- 
ping disorder sd, T — 0.05f/ and an on-site speckle disorder 
strength A = U. Main 3D figure: disorder-averaged MFP tl) 
for speckle on-site and experimentally corresponding hopping 
disorder. For orientation purposes, the MF phase boundary 
for the pure system is plotted as a thin dashed white line. 
The value J/U refers to the most likely value of the hopping 
disorder distribution pj{J), whereas the width of the distri- 
bution Pj{J) is constant. Upper left inset: comparison of the 
order parameter at fixed jj, — 2.6U with and without hopping 
disorder, clearly demonstrating that the addition of hopping 
disorder stabilizes the SF phase and leads to a lower criti- 
cal value of J in this regime. Upper right inset: Subsequent 
comparison of phase borders in the n/U-J/U-plane. 



^tot = ^/wi + wi^g. (31) 

For every fixed value of the lattice height s, J and U 
are extracted from the single particle Wannier function 
(beyond the approximation for deep lattices^), the to- 
tal trapping frequency (|31]) is calculated and the local 
effective trap energies are assigned on a sufficiently large 
3D lattice. The chemical potential /x is adjusted using 
Ridder's method^ to obtain a specified total particle 
number. The condensate fraction, depicted in Fig. [T4l 
is subsequently averaged over the local values obtained 
by SMFT, weighted by the respective density. 



After ramping the lattice to a sufficiently high intensity, 
the entropy can be calculate up to first order (neglecting 
terms containing J)^ 



1 

N 



li-i{E{M)) + I3E 



(28) 



where M is the number of sites and S(M) is the 
grand canonical partition function. Equating these two 
expressions for a sufficiently high initial temperature 
^sTinitiai > O.Q5Eii leads to the rclatior 



knTt 



B-t final 



U 

3^ 



(fesT..,.,.,, + 0.177^flJ 



(29) 



For *^Rb in an optical with a wave length of A = 812nm 
and intensity s = llEfi, the disorder-averaged interac- 
tion constant is U/Er = 0.355. As a typical, conserva- 
tive estimate of the initial temperature before the lat- 
tice ramp-up in current experiments^ we use the value 
^initial = 0.13/iX, for which the relation (|29p predicts a 
final temperature of kgT w O.llEjj = 0.32U after the 
ramp-up. 
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FIG. 14: (Color online) The condensate fraction as a function 
of the lattice strength s for various total particle numbers in a 
harmonic trap, calculated within LDA and SMFT (connected 
dotted lines). The calculations were performed at fixed tem- 
perature kT = 0.3(7 and disorder strength A = If/. 

We used the following values for the experimental pa- 
rameters: wq ~ 110 fim, the lattice laser wavelength was 
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set to Aiat = 812nm, and the magnetic trapping frequency 
a;„,.,^ = 2'k40Hz. 

In contrast to the behavior of the order parameter at 
the transition point, the condensate fraction does not 
follow a power law decay in the finite trap, but is smeared 
out a the transition point, as superfluid regions in the 
trap decay in size with increasing s. 



VII. CONCLUSION 

We have described the stochastic mean-field theory in 
detail on a methodological level, and extend it to incor- 
porate finite temperature effects. Subsequently, we have 
applied it to ultracold atoms in an optical lattice with 
uncorrclated on-site box disorder distribution and dis- 
cussed the intricate interplay between interaction, tun- 
neling energy, disorder, filling and finite temperature ef- 
fects. Furthermore, we have presented, to the best of 
our knowledge, the first quantitative theoretical calcula- 
tions for speckle disorder, which leads to a qualitatively 
different phase diagram than for box disorder and are 
of immediate experimental relevance^^. For this case, 
we have discussed the characteristic features of the vari- 
ous phases and presented phase diagrams, both at fixed 
chemical potential and at fixed density. Below a critical 
temperature, we find disorder-induced condensation and 
multiple reentrant behavior, both for box and speckle 
disorder. The temperature in recent experimentsi^ is es- 
timated and found to be too high yet to observe disorder- 
induced condensation. We also find that including hop- 
ping disorder in addition to local on-site disorder for a 
realistic distribution of speckle parameters enhances the 
insulator and jumps in the order parameter within the SF 
phase indicate a series of transitions. An LDA-I-SMFT 
calculation has been performed to incorporate the effects 
of an external trap for on-site speckle disorder. 
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Appendix A: Multi-Grid discretization 



behavior in the vicinity of the SF-insulator transition, as 
well as simultaneously correctly describing the distribu- 
tion at large values {tp ^ 1). To circumvent this problem, 
we use a superposition of equidistant grids, which enables 
us to use the FFT algorithm on each of these grids indi- 
vidually (see Fig. [T5|). This procedure relies on the fol- 
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final "f/^-grid 



FIG. 15; Illustration of the multi-grid procedure. 

lowing property of the convolution the Z-fold convo- 
lution of the truncated function Ptiip) — ©(a — "0) Pii') 
with P{ijj) ~ for -0 < and a < ip„^^^/Z is identical 
to the Z-fold convolution of P{ip) up to an easily de- 
terminable normalization constant on the interval [0, a]. 
Here '(/'max is the largest value of the grid, if this is fi- 
nite (as in the discretized numerical case). In our cal- 
culations the number of grids used (typically 1 ... 6) 
with 200 . . . 1000 points per grid is adjusted dynamically 
within the iteration procedure, depending on the posi- 
tion of the most likely value of P{ip) and the convergence 
properties. 



Appendix B: Numerical calculation of the CDF 

F{ip\ri) for box disorder 

We found the numerically most efficient method for 
calculating and tabulating the two-dimensional CDF 
F{ip\ri) to be the following: 

For every fixed value of rj, consider the function ([7]) 
g{fj,','r]) on the interval /i' = (/i — e) £ [/i — A/2,/i + 
A/2]. We define the local minima and maxima of this 
function (including the end points) in increasing order as 
{^Lt"'", /^^'", . . .} and {a*""'', /Lt^"", . . .}. Furthermore the 
n-th monotonically increasing and decreasing function on 
the restricted interval, provided that this interval exists, 
is denoted by 

g("-")(A.',77)=5(/i',^) for /i' € (a^^, Ar") (Bl) 
g(--")(A.',77)=5(/i',^) for M'e(A^^^A;^) (B2) 



where 



f^r'^ = min{pi»="'|/x;';;"" > /x™} (B3) 

rn 

f^n'" = min{Ai™|Mm" > ^r"} (B4) 



To employ the FFT algorithm an equidistant grid is re- 
quired, which is not compatible with having a very high 
resolution at small values of t/j (^ 10~^) to capture the 



By construction, these functions are invertible in /i' 
within the defined range, allowing us to introduce the 
functions 
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if ■0 < gifJ-n^iV) 
if e {g{fi';^'",ri), 

if > gi.K'^^-.v) 

if V' < 9{fj-n"\V) 

if ■(/; e {g{fl';^'",ri), g(Air">^)) 

if ■iA > g(/^r":'7) 



(B5) 
(B6) 



in terms of which the CDF can be written as a super- 
position 



Appendix C: PDF for a product of random variables 

Inchiding hopping disorder into the SMFT leads to the 
task of having to calculate the probability distribution for 
the newly defined random variable (j) = •^V'l if the PDFs 
for J and ijj are known functions pj{ J) and This 
can be done by taking the logarithm s = ln(0) = x + y 
and using the convolution theorem for the PDFs {x) ~ 

Pj{e^) (analogously for Py{y)) of the random variables 
X = ln(J) and y = ln(?/'). The resulting distribution for 
(j) can be denoted in compact form as 



P^(^) = J dxPj{e^)P{ 



-1- 



(CI) 



Using the FFT-algorithm for the numerical computation 
of this operation leads to a vast increase in performance, 
however it requires the functions Px{x) and Py{y) to be 
interpolated on an equidistant grid. 



Appendix D: Compressibility in the insulating 
phases 

In the special case with ^("0) = S{ip), the disorder in- 
tegration can be performed explicitly and expressions for 
the density and the compressibility can be found. Using 
partial integration and the property limg^oo p(ie) — 
for the on-site probability distribution, the disorder av- 
eraged density at finite temperature can be expressed as 



niP,U,^l,A)^ / dep{e) 



m=0 ' 



-/3£;™(C/,M-e) 



Z_^m— 



\ 771 

\ m 



1 f^^dpie 



P J de 



(Dl) 



with E„i{U, /i') = U mirn — l)/2 — /i'm. 
In the specific case of box disorder with 

^ = i[<5(6 + A/2)-<5(e-A/2)] 
the disorder averaged density takes on the form 



n(/3, U, A) = — In 



g-/3£;„.(M-A/2) 



(D2) 



(D3) 



and the compressibility = ^ can be directly evaluated 



1 fEr 



m e 



-/3iS„(p-A/2) 



-/3£;™(ai+A/2) 



(D4) 



For speckle disorder, on the other hand, one has 

9p(e) ^ ^ _ e(e) _,/^ 
de A A2 ^ 

and the compressibility takes on the form 



(D5) 



1 

A _ 

^ />oo 



m e 



-l3Em{fl) 







(D6) 



Appendix E: Local Green's functions and DOS 

To obtain the local DOS in an insulating phase, we 
calculate the single particle Green's functions 



G>{t)^mh\o)) 
G<{t)^{bHmt)) 

Git) = -^[e(^) G>{t) + Qi-t) G<(i)] 



(El) 



for local Fock states at finite temperature, where b{t) is 
the on-site particle annihilation operator in the Heisen- 
berg representation. The Fourier transformed Green's 
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function can be calculated and takes on the form 
G{uj) = / dte"^*G{t) 

(■^m(m — 1) — /j'm) 
UJ — Um + ll' ~ i"f 



= lim } e 

^' ' 771—0 



(E2) 



— C/(m — 1) + /z' — 27 



(w - C/(m - 1) + //')^ + 7^ 



where Z{fi) = ^^^^ e-/3(f m(m-i)-M™) jg the local par- 
tition function. This is related to the single particle DOS 
by 



1_ 



= — Im(GM) 



^ OO 

x6{iu - Um + /i') + TO(5(cj - C/(m - 1) + ^f')] • 

(E3) 

Averaging over the on-site energy distribution for 
speckle disorder leads to the final expression 

p{uj,fi,A,P)= / dep{e) p{uj,fi~ e) 



_1 ^ (m + l) e(a; - Um + p) 
~ A ^ Z{Um - uj) 

m=0 ^ ' 

X e~ a""^" -/3(m(^-f m(m+l)) ^ ^ 

(E4) 

whereas for a box disorder distribution one obtains^ 

p(w,/X,A,/3) = — 2^(TO+l)- 



Jn=0 



Z{Um — uj) 



X e 



(E5) 



Appendix F: Method of incorporating thermal 
fluctuations explicitly into SMFT 

Instead of performing the thermal average before con- 
structing the conditional probability distribution for ip, 
the SMFT furthermore also allows to explicitly facilitate 
the thermal fluctuations of ip in the probability distribu- 
tion. We have not yet performed the numerical calcula- 
tion, but will outline the procedure. 



Pii'lv) = ^ / dep{e) e 



Tr(6( 



Tr (e-'^"<-^-^)) 



, (Fl) 



where &{x) is the Heaviside step function. 

This can be formulated in the following way: for 
fixed external parameter rj (hence conditional probabil- 
ity distribution) the on-site energy is randomly drawn 
from p{e) and the resulting single site Hamiltonian 
is diagonalized. For the different eigenstates \i{e,T])) 
with eigenenergies Ei the respective expectation value 
{i{e,r])\ b\i{e,r])) is calculated. Within the grand canon- 
ical ensemble, this expectation value has the probability 
j-rpj. (•g-/3//(,,,e)^j-i g-/3E. of occurring. This probability 
is uncorrelated to the probability of a certain on-site en- 
ergy e occurring. Expressing this mathematically, we ob- 
tain the explicit formula for the conditional probability 
distribution 



Tr e-^"^^''~> 



Pi^lv) =^ / depie) 



(F2) 



Equivalently on a formal level, one may work with the 
cumulative conditional probability function 



F{ij\r]) = I dep{e) Tr (^e''^"''''-''^^ 



(F3) 



with 



g,{p-e,-q) = {i{e,ri)\b\i{e,if)) 
To evaluate the expression 



(F4) 



de/,(e)e(V'-<?.(/^ -£,»?)) (F5) 



we perform a variable substitution 

e ^ Xi{e) 

and subsequently 

dx,{e) 
dxi = — ; — de. 



(F6) 



(F7) 



Using the central theorem of calculus one can explicitly 
construct 



Since 



Xi{e) 



de'Me'). 



g-/3B,(e,77) 

Tr(e-^^(^.'')) 



(F8) 



(F9) 
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is a non-ncgativc function the monotonously increasing 
function Xi{e) O ei{x) is invertible in the relevant range. 
The conditional cumulative density function can then be 
explicitly calculated from 

Fii'lv) = ^ / 9(7^ - g^ifi ~ e,(x)), r,) (FIO) 

2 — 1 

As in the SMFT for T = 0, the self-consistency condition 
reads 

m = d7jQ{ij)F{4,\Tj), (Fll) 

dyj Jo 



where Q{r]) is the fold convolved and rescaled function 
of P{ip) with the random variable r] = J'Yld=i'4'i- The 
conditional cumulative density F{il)\vi) now also contains 
the thermal fluctuations explicitly in the distribution (i.e. 
they are not averaged over within the self-consistency 
loop), so the subsequent determination of P{tp) is iden- 
tical to the previous cases. 
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